Tanzania
Rasmus Benestad
May 10, 2016
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
## Script to downscale temperature from Tanzania.
## Data from Habiba.
library(esd)Functions
Prepare the predictands
predictands <- function(param='tmin',incl.ghcn=FALSE) {
## Prepare the predictand data
#print("predictands")
fname <- paste(param,'.tz.rda',sep='')
path <- paste('Habiba/',param,sep='')
if (!file.exists(fname)) {
files <- list.files(path,pattern='rda',full.names=TRUE)
locs <- list.files(path,pattern='rda')
dse <- grep('dse.tm',files)
files <- files[-dse]; locs <- locs[-dse]
print(files)
locs <- substr(locs,1,nchar(locs)-4)
Y <- NULL
for (i in 1:length(files)) {
load(files[i])
y1 <- eval(parse(text=locs[i]))
if (is.null(Y)) Y <- y1 else
Y <- combine.stations(Y,y1)
}
save(file=fname,Y)
plot(Y)
} else load(fname)
#print('GHCN')
if (incl.ghcn) {
xstations <- switch(param,'tmin'='Habiba/tn.eafrica.rda',
'tmax'='Habiba/tx.eafrica.rda')
load(xstations)
if (param=='tmin') y <- tn.eafrica else y <- tx.eafrica
y <- subset(y,it=c(1950,2015))
Y <- combine(Y,y)
}
map(Y,FUN='mean',cex=-2)
invisible(Y)
}Function that applies the downscaling to a combination of predictands in terms of principal component analysis (PCA) results and predictors that involve large multi-model ensembles.
downscale <- function(Y,predictor,it='djf',param='t2m',FUN='mean',FUNX='mean',
period=c(1950,2015),plot=FALSE,rcp='rcp45',verbose=FALSE,
lon=c(0,100),lat=c(65,90),eofs=1:6,n=10,rel.cord=FALSE,select=NULL) {
print('downscale')
## Use a time and space window:
Y <- subset(Y,it=period)
Y <- subset(Y,is=list(lat=lat,lon=lon))
## Estimate seasonal means & weed out stations with little data
Y4 <- subset(as.4seasons(Y,FUN=FUN,nmin=60),it=it)
ok <- apply(coredata(Y4),1,nv)
Y4 <- subset(Y4,it=ok>0)
nok <- apply(coredata(Y4),2,nv)
Y4 <- subset(Y4,is=nok>15)
print(paste(round(100*sum(!is.finite(Y4))/length(!is.finite(Y4))),'% missing',sep=''))
if (plot) map(Y,FUN=FUN,cex=-2)
nmiss <- round(100*sum(!is.finite(Y4))/length(!is.finite(Y4)))
print(paste(nmiss,'% missing',sep=''))
## Fill missing data using PCA-based regression
Z <- pcafill(Y4)
pca <- PCA(Z,n=n)
if (plot) plot(pca)
## Downscale results
print('DSensemble')
dse.pca <- DSensemble(pca,predictor=predictor,FUNX=FUNX,verbose=verbose,
biascorrect=TRUE,rcp=rcp,eofs=eofs,select=select,
lon=lon,lat=lat,rel.cord=rel.cord,it=it)
attr(dse.pca,'N.missing') <- nmiss
invisible(dse.pca)
}Function that extracts results for a given year
distill <- function(x,it=2050) {
y <- subset(x,it=it)
attributes(y) <- NULL
y <- mean(y)
return(y)
}Function that grids the results
gridmap <- function(dse.results,param=NULL,
fast=TRUE,it=c(2070,2099),it0=c(1979,2012),verbose=FALSE,
breaks=seq(0,7,length.out=15),pal='warm',rev=FALSE) {
if (verbose) print('gridmap')
library(LatticeKrig)
if (is.null(param)) param <- varid(dse.results)
if (verbose) print ('as.station')
if (inherits(dse.results,'pca')) dse.station <- as.station(dse.results) else
dse.station <- dse.results
if (verbose) print('distill')
Z0 <- lapply(dse.station,distill,it=it0)
Z <- lapply(dse.station,distill,it=it)
lon <- unlist(lapply(dse.station,lon))
lat <- unlist(lapply(dse.station,lat))
alt <- unlist(lapply(dse.station,function(x) alt(attr(x,'station'))))
X <- attr(dse.station[[1]],'station')
z0 <- matrix(unlist(Z0),1,length(dse.station))
z <- matrix(unlist(Z),1,length(dse.station))
data(etopo5)
etopo5 <- subset(etopo5,is=list(lon=range(lon),lat=range(lat)))
etopo5[etopo5<=0] <- NA
if (verbose) print('grid')
## Set the grid to be the same as that of etopo5:
grid <- structure(list(x=lon(etopo5),y=lat(etopo5)),class='gridList')
#3 Flag dubplicated stations:
ok <- !(duplicated(lon) & duplicated(lat))
## Spread in the 90-percente interval changing
obj <- LatticeKrig( x=cbind(lon[ok],lat[ok]),
y=z[1,ok] - z0[1,ok],Z=alt[ok] )
## obj <- LatticeKrig( x=cbind(lon[ok],lat[ok]), y=z[2,ok],Z=alt[ok])
w <- predictSurface(obj, grid.list = grid,Z=etopo5)
w$z[is.na(etopo5)] <- NA
## Get rid of packages that have functions of same name:
detach("package:LatticeKrig")
detach("package:fields")
detach("package:spam")
detach("package:grid")
detach("package:maps")
## Convert the results from LatticeKrig to esd:
W <- w$z
attr(W,'variable') <- param
attr(W,'unit') <- 'degC'
attr(W,'longitude') <- w$x
attr(W,'latitude') <- w$y
class(W) <- class(etopo5)
## Make a projection that zooms in on the Barents region
#dev.new()
colbar <- list(breaks=round(breaks,2),rev=rev,pal=pal)
if (verbose) print(colbar)
if (is.null(rev)) rev <- switch(param,'t2m'=FALSE,'pr'=TRUE)
Wx <- max(abs(W),na.rm=TRUE)
if (is.null(breaks)) breaks <- round(seq(0,Wx,length=31),2)
if (fast)
map(W,xlim=range(lon(W)),ylim=range(lat(W)),
colbar=colbar,new=FALSE)
else {
attr(W,'variable') <- NULL
attr(W,'unit') <- NULL
map(W,xlim=range(lon(W)),ylim=range(lat(W)),projection='sphere',
colbar=colbar,verbose=verbose)
}
figlab(paste(it,collapse='-'),ypos=0.99)
}Function that applies tests to theresults through applying PCA
pcaasses <- function(Y) {
## Cut away blocks with missing data:
nv <- apply(coredata(Y),2,nv)
Y <- subset(Y,is=nv > 12000)
Y4s <- as.4seasons(Y,nmin=60)
y1 <- subset(Y4s,it='Jan')
y2 <- subset(Y4s,it='Apr')
y3 <- subset(Y4s,it='Jul')
y4 <- subset(Y4s,it='Oct')
## Fill in missing
y1 <- pcafill(y1)
y2 <- pcafill(y2)
y3 <- pcafill(y3)
y4 <- pcafill(y4)
## Apply PCA:
pca.djf <- PCA(y1)
pca.mam <- PCA(y2)
pca.jja <- PCA(y3)
pca.son <- PCA(y4)
## Examine the results:
print('DJF - PCA')
plot(pca.djf,new=FALSE)
print('MAM - PCA')
plot(pca.mam,new=FALSE)
print('JJA - PCA')
plot(pca.jja,new=FALSE)
print('SON - PCA')
plot(pca.son,new=FALSE)
## compare the station data with reanalysis - based on CCA
## Read reanalysis:
X <- retrieve('air.mon.mean.nc',lon=range(lon(Y))+c(-3,3),lat=range(lat(Y))+c(-3,3))
X4s <- as.4seasons(X)
eof.djf <- EOF(subset(X4s,it='Jan'))
eof.mam <- EOF(subset(X4s,it='Apr'))
eof.jja <- EOF(subset(X4s,it='Jul'))
eof.son <- EOF(subset(X4s,it='Oct'))
cca.djf <- CCA(pca.djf,eof.djf)
cca.mam <- CCA(pca.mam,eof.mam)
cca.jja <- CCA(pca.jja,eof.jja)
cca.son <- CCA(pca.son,eof.son)
print('DJF - CCA')
plot(cca.djf,new=FALSE)
print('MAM - CCA')
plot(cca.mam,new=FALSE)
print('JJA - CCA')
plot(cca.jja,new=FALSE)
print('SON - CCA')
plot(cca.son,new=FALSE)
## The quick and effiecint way to downscale all the stations:
ds.djf <- DS(subset(pca.djf,pattern=1:4),eof.djf)
ds.mam <- DS(subset(pca.mam,pattern=1:4),eof.mam)
ds.jja <- DS(subset(pca.jja,pattern=1:4),eof.jja)
ds.son <- DS(subset(pca.son,pattern=1:4),eof.son)
print('DJF - DS')
plot(ds.djf,new=FALSE)
#dev.new()
print('MAM - DS')
plot(ds.mam,new=FALSE)
#dev.new()
print('JJA - DS')
plot(ds.jja,new=FALSE)
#dev.new()
print('SON - DS')
plot(ds.son,new=FALSE)
return(list(pca.djf=pca.djf,pca.mam=pca.mam,pca.jja=pca.jja,
pca.son=pca.son,y.djf=y1,y.mam=y2,y.jja=y3,y.son=y4))
}Maximum temperature
The code below shows how the functions are used and applied to the data for the daily maximum temperature \(T_x\). The first chunck of code carries out some preliminary analysis, the second shows how the results are derived, and the third show how the plots are made.
#-----------------------------------------------------------------------
# Define season and parameter
param <- 'tmax'
FUN <- 'mean'
#reanalysis <- 'air.mon.mean.nc'
reanalysis <- 'ERAINT_t2m_mon.nc'
FUNX <- 'mean'
data.path <- 'dse.Tz.ERAINT/'
predictands(param) -> Y
Y <- subset(Y,it=c(1961,2011))
lon=range(lon(Y))+c(-10,10); lat=range(lat(Y))+c(-15,15)
## Check the predictands:
X <- pcaasses(Y)## [1] "DJF - PCA"
## [1] "MAM - PCA"
## [1] "JJA - PCA"
## [1] "SON - PCA"
## [1] "Warning : Calendar attribute has not been found in the meta data and will be set automatically."
## [1] "DJF - CCA"
## [1] "MAM - CCA"
## [1] "JJA - CCA"
## [1] "SON - CCA"
##
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%[1] "DJF - DS"
## [1] "MAM - DS"
## [1] "JJA - DS"
## [1] "SON - DS"
The PCA suggest that most of the variance in the station data for daily maximum temperature \(T_x\) can be represented by a small number of leading modes. For the SON season, however, the second mode is more prominent, suggesting more complicated structure where spatial variability has less coherency between the stations.
The CCA suggests dissimilar spatial temperature pattern with the higest temporal correlation in the predictand and predictor. The closest spatial similarity is seen in SON.
The DS analysis applied to the reanalysis suggests strongest match in MAM and weakest in DJF. The spatial patterns differ to some extent and there is not a great match between the predictor and predictand.
## Get the predictand -> Y
for (it in c('djf','mam','jja','son')) {
## Get the large-scale predictor:
if (!exists('predictor')) {
T2M <- retrieve(reanalysis,lon=lon,lat=lat)
predictor <- subset(as.4seasons(T2M,FUNX=FUNX),it=it)
} else if (length(month(subset(predictor,it=it)))==0)
predictor <- subset(as.4seasons(T2M,FUNX=FUNX),it=it)
#print(paste('Generating dse.',param,'.tz.rcp45.',it,'.rda',sep=''))
## Carry out the downscaling:
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp45.',it,'.rda',sep=''))) {
#dev.new()
dse.t2m.tz.rcp45 <- downscale(Y,predictor,it,param,FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp45.',it,'.rda',sep=''),dse.t2m.tz.rcp45)
}
#print(paste('Generating dse.',param,'.tz.rcp26.',it,'.rda',sep=''))
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp26.',it,'.rda',sep=''))) {
dse.t2m.tz.rcp26 <- downscale(Y,predictor,it,param,rcp='rcp26',
FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp26.',it,'.rda',sep=''),dse.t2m.tz.rcp26)
}
#print(paste('Generating dse.',param,'.tz.rcp85.',it,'.rda',sep=''))
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp85.',it,'.rda',sep=''))) {
dse.t2m.tz.rcp85 <- downscale(Y,predictor,it,param,rcp='rcp85',
FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp85.',it,'.rda',sep=''),dse.t2m.tz.rcp85)
}
}
#print('--- Completed downscaling ---')The March-May season.
## analysis:
load(paste(data.path,'dse.',param,'.tz.rcp85.mam.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
plot(subset(z,is='Arusha'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('figure4g.',param,'.pdf'))
plot(subset(z,is='Musoma'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('figure4h.',param,'.pdf'))
gridmap(z)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.pdf',sep=''))
validate(z)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.diagnose.pdf',sep=''))The plume plots look OK for Arusha and Musoma, and the observations fit reasonably well within the envelope of downscaled ensemble (RCP8.5). The gridded map suggests strongest warming in the north (~5C increase in 2099).
The rank test (Wilcox) indicates reasonable scores, although the rank distribution is not entirely flat. The diagnostic (target plot) for the diagnosis suggest a reasonable reproduction of trend (horizontal axis) and variance (vertical axis), although four stations were associated with higher variance than suggested by the ensemble and the observed trends were weaker than suggested by the ensemble.
The December-February season
load(paste(data.path,'dse.',param,'.tz.rcp85.djf.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.pdf',sep=''))
validate(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.validate.pdf',sep=''))
diagnose(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.diagnose.pdf',sep=''))The projected warming was again hiher in the north (~3C), and the evaluation suggested reasonable skill, although not perfect.
The June-August season
load(paste(data.path,'dse.',param,'.tz.rcp85.jja.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.pdf',sep=''))
validate(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.validate.pdf',sep=''))
diagnose(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.diagnose.pdf',sep=''))In June-August, the projected warming in \(T_x\) was found to be greatest the west (~3C), but the assessment of the rank score suggested a poor match between observations and the downscaled ensemble.
The September-November season
load(paste(data.path,'dse.',param,'.tz.rcp85.son.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.pdf',sep=''))
validate(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.validate.pdf',sep=''))
diagnose(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.diagnose.pdf',sep=''))The results suggest strongest warming in the north, but the map is unrealistic with only a latitude-based gradient. The assessment of the ranked score also suggests a weak performance.
Test of PCA
Check a random station: is it well-represented by the PCA?
## Check the PCA-results:
z.djf <- as.station(X$pca.djf)
y.djf <- X$y.djf
par(new=FALSE)
plot.zoo(subset(z.djf,is=1),subset(y.djf,is=1),
main=paste('Test: Winter temperature at',loc(subset(y.djf,is=1))))
plot(subset(z.djf,is=1),lwd=3,col='grey70')
lines(subset(y.djf,is=1))The data for Arusha was recovered from the PCA products.
Minimum temperature
Get and prepare the predictand data for the daily minimum temperature \(T_n\).
#-----------------------------------------------------------------------
# Define season and parameter
param <- 'tmin'
predictands(param) -> Y
Y <- subset(Y,it=c(1961,2011))
lon=range(lon(Y))+c(-10,10); lat=range(lat(Y))+c(-15,15)
## Check the predictands:
X <- pcaasses(Y)## [1] "DJF - PCA"
## [1] "MAM - PCA"
## [1] "JJA - PCA"
## [1] "SON - PCA"
## [1] "Warning : Calendar attribute has not been found in the meta data and will be set automatically."
## [1] "DJF - CCA"
## [1] "MAM - CCA"
## [1] "JJA - CCA"
## [1] "SON - CCA"
##
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%[1] "DJF - DS"
## [1] "MAM - DS"
## [1] "JJA - DS"
## [1] "SON - DS"
The PCA suggests that most of the variance in \(T_n\) was captured by a small number of leading modes. The CCA found that the patterns in the predictand and predictor with highest correlation were not so similar in terms of spatial structure.
The DS applied to the reanalysis found highest scores for MAM, and again, the spatial patterns were dissimilar.
Apply the ESD models to implement the downscaling
## Get the predictand -> Y
for (it in c('djf','mam','jja','son')) {
## Get the large-scale predictor:
if (!exists('predictor')) {
T2M <- retrieve(reanalysis,lon=lon,lat=lat)
predictor <- subset(as.4seasons(T2M,FUNX=FUNX),it=it)
} else if (length(month(subset(predictor,it=it)))==0)
predictor <- subset(as.4seasons(T2M,FUNX=FUNX),it=it)
print(paste('Generating dse.',param,'.tz.rcp45.',it,'.rda',sep=''))
## Carry out the downscaling:
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp45.',it,'.rda',sep=''))) {
#dev.new()
dse.t2m.tz.rcp45 <- downscale(Y,predictor,it,param,FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp45.',it,'.rda',sep=''),dse.t2m.tz.rcp45)
}
print(paste('Generating dse.',param,'.tz.rcp26.',it,'.rda',sep=''))
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp26.',it,'.rda',sep=''))) {
dse.t2m.tz.rcp26 <- downscale(Y,predictor,it,param,rcp='rcp26',
FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp26.',it,'.rda',sep=''),dse.t2m.tz.rcp26)
}
print(paste('Generating dse.',param,'.tz.rcp85.',it,'.rda',sep=''))
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp85.',it,'.rda',sep=''))) {
dse.t2m.tz.rcp85 <- downscale(Y,predictor,it,param,rcp='rcp85',
FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp85.',it,'.rda',sep=''),dse.t2m.tz.rcp85)
}
}## [1] "Generating dse.tmin.tz.rcp45.djf.rda"
## [1] "Generating dse.tmin.tz.rcp26.djf.rda"
## [1] "Generating dse.tmin.tz.rcp85.djf.rda"
## [1] "Generating dse.tmin.tz.rcp45.mam.rda"
## [1] "Generating dse.tmin.tz.rcp26.mam.rda"
## [1] "Generating dse.tmin.tz.rcp85.mam.rda"
## [1] "Generating dse.tmin.tz.rcp45.jja.rda"
## [1] "Generating dse.tmin.tz.rcp26.jja.rda"
## [1] "Generating dse.tmin.tz.rcp85.jja.rda"
## [1] "Generating dse.tmin.tz.rcp45.son.rda"
## [1] "Generating dse.tmin.tz.rcp26.son.rda"
## [1] "Generating dse.tmin.tz.rcp85.son.rda"print('--- Completed downscaling ---')## [1] "--- Completed downscaling ---"The March-May season
## analysis:
load(paste(data.path,'dse.',param,'.tz.rcp85.mam.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
plot(subset(z,is='Arusha'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('figure4g.',param,'.pdf'))
plot(subset(z,is='Musoma'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('figure4h.',param,'.pdf'))
gridmap(z)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.pdf',sep=''))
validate(z)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.diagnose.pdf',sep=''))The plume plots for Arusha and Musoma look ok. Strongest warming (~3C) is projected for the eastern part of the country. THe assessment of the downscaled results in terms of ranked score suggest that the observations have a different character to the downscaled results (cold bias).
The December-February season
load(paste(data.path,'dse.',param,'.tz.rcp85.djf.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.pdf',sep=''))
validate(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.diagnose.pdf',sep=''))The projections suggest strongest warming in the southeast (~2C), and again, the evaluation indicates a mismatch between observations and downscaled results. Poor performance.
The June-August season
load(paste(data.path,'dse.',param,'.tz.rcp85.jja.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.pdf',sep=''))
validate(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.diagnose.pdf',sep=''))Strongest warming in the southwest (~3C), and higher ranked scores than the two previous seasons, but the observations have a higher variance than accounted for by the downscaled ensemble and the the trends in the downscaled results are weaker than those observed.
The September-November season
load(paste(data.path,'dse.',param,'.tz.rcp85.son.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.pdf',sep=''))
validate(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.diagnose.pdf',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
plot(subset(z,is='Arusha'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('SON')
plot(subset(z,is='Musoma'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('SON')
Greatest warming in the southwest, but these results look suspect. The plume-plot for Arusha exhibits a spread that increases over time that spans from cooling to warming.
Check a random station: is it well-represented by the PCA?
## Check the PCA-results:
z.djf <- as.station(X$pca.djf)
y.djf <- X$y.djf
plot.zoo(subset(z.djf,is=1),subset(y.djf,is=1),
main=paste('Test: Winter temperature at',loc(subset(y.djf,is=1))))
plot(subset(z.djf,is=1),lwd=3,col='grey70')
lines(subset(y.djf,is=1))The original data is recovered from the PCA product.
Summary
There are problems with some stations, which is visible in the PCA, CCA and DS carried out with just reanalysis. The spatial patterns indicate non-smooth and incoherent structures, suggesting problems with some of the stations and a mismatch with the reanalysis. Furthermore, the cross-validation reveal moderate scores (r typically less than 0.7), and the evaluations based on trends, range of inter-annual variability and rank score reveal some discrepancies.
The ERAINT is probably too short for proper model calibration.